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ABSTRACT 


Our  research  program  consists  of  four  components,  each  involving  some  aspect  of  multiple-event  analysis:  (1)  high- 
precision  waveform  cross-correlation  (WCC)  for  arrival  time  estimation,  (2)  robust  event  clustering,  (3)  waveform 
decomposition  and  source  wavelet  deconvolution  reshaping,  (4)  double-difference  (DD)  multiple-event  location  and 
tomography.  Our  research  focused  initially  on  the  development  and  testing  of  these  seismic  analysis  methods  using 
"ground-truth"  (GT)  datasets  at  different  scales  (local  and  regional),  and  is  being  followed  by  the  application  of 
these  methods  to  a  "test-bed"  regional-distance  seismic  dataset,  including  tests  on  real-time  and  simulated  real-time 
data  streams.  Here  we  highlight  some  of  the  more  significant  project  achievements  from  components  1, 2,  and  4. 

(1)  Differential  arrival  times  for  pairs  of  seismic  events  observed  at  the  same  station  are  often  calculated  by  WCC 
techniques.  Researchers  often  choose  the  differential  times  to  use  based  on  the  associated  cross-correlation  (CC) 
values  exceeding  a  specified  threshold.  When  two  similar  time  series  are  corrupted  by  correlated  noise,  the  time 
delay  estimate  calculated  with  the  WCC  technique  may  not  be  reliable.  Thus  the  selection  of  a  threshold  value  is 
important.  If  it  is  set  too  high,  then  only  a  limited  number  of  very  accurate  differential  time  data  are  available  to 
constrain  the  relative  positions  of  earthquakes.  If  the  threshold  value  is  set  too  low,  then  many  unreliable  differential 
time  estimates  are  used  and  they  will  negatively  affect  the  final  relocation  results.  The  bispectrum  method  can 
suppress  the  correlated  Gaussian  noise  sources  in  two  similar  time  series  and  be  used  to  obtain  the  relative  time 
delay  between  them.  We  compute  time  delay  estimates  between  the  waveform  pairs  with  the  bispectrum  method  and 
use  them  to  verify  the  WCC-determined  one.  This  technique  can  provide  quality  control  over  the  selected  time  delay 
estimates  and  potentially  provide  more  differential  travel  time  data  for  close  events  pairs,  by  verifying  the  reliability 
of  differential  times  that  might  not  meet  the  threshold  criterion. 

(2)  We  have  tested  a  two-step  event  location  technique  with  a  California  GT  data  spt.  In  step  1,  event  locations  are 
estimated  via  WCC  of  Hilbert-envelope  waveforms.  The  Step  1  procedure  accurately  located  the  epicenters  within 
20  km  (approximately  1  bin)  of  the  NCSN  catalog  location  for  16  of  the  20  test  events.  However,  depth  estimates 
were  poorly  constrained.  The  Step  2  refinement  based  on  WCC  of  the  high-frequency  P  arrivals  followed  by  DD 
relocation  improved  both  epicenter  and  depth  estimates.  This  two-step  procedure  offers  advantages  over 
conventional  location  techniques  including  a  Step  1  event  location  not  dependent  on  the  accuracy  of  catalog  picks 
and  a  Step  2  improvement  in  location  accuracy  using  the  high-accuracy  DD  technique. 

(3)  We  have  developed  the  new  method  of  DD  tomography.  The  inversion  code  tomoDD  takes  the  double¬ 
difference  location  method  (hypoDD)  a  step  further  to  solve  for  3D  velocity  structure  simultaneously  with 
hypocenter  locations  using  both  catalog  picks  and  differential  arrival  times  (WCC  and  catalog).  Compared  to 
conventional  tomography,  the  result  is  a  sharpening  of  the  seismicity  distribution  equal  in  quality  to  hypoDD  plus  a 
sharpening  of  the  velocity  image  due  to  removal  of  the  majority  of  the  location  scatter.  We  have  developed  local- 
and  regional-scale  versions  of  tomoDD,  and  have  applied  each  to  synthetic  and  real  datasets. 
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OBJECTIVE 


Our  research  program  consists  of  four  components,  each  involving  some  aspect  of  multiple-event  analysis:  (1)  high- 
precision  waveform  cross-correlation  (WCC)  for  arrival  time  estimation,  (2)  robust  event  clustering,  (3)  waveform 
decomposition  and  source  wavelet  deconvolution  reshaping,  (4)  double-difference  (DD)  multiple-event  location  and 
tomography.  Our  research  focused  initially  on  the  development  and  testing  of  these  seismic  analysis  methods  using 
"ground-truth"  (GT)  datasets  at  different  scales  (local  and  regional),  and  is  being  followed  by  the  application  of 
these  methods  to  a  "test-bed"  regional-distance  seismic  dataset,  including  tests  on  real-time  and  simulated  real-time 
data  streams.  Here  we  highlight  some  of  the  more  significant  project  achievements  from  components  1,  2,  and  4. 

RESEARCH  ACCOMPLISHED 


Bispectrum  analysis 

Differential  arrival  times  for  pairs  of  seismic  events  observed  at  the  same  station  are  often  calculated  by  WCC 
techniques.  These  differential  times  (or  adjusted  picks  derived  from  them)  can  be  used  to  improve  dramatically  the 
earthquake  relocation  results  (Got  et  al.,  1994;  Dodge  et  ah,  1995;  Shearer,  1997;  Rubin  et  ah,  1999;  Waldhauser 
and  Ellsworth,  2000;  Schaff  et  ah,  2002).  Researchers  often  choose  the  differential  times  to  use  based  on  the 
associated  CC  values  exceeding  a  specified  threshold.  For  example,  Schaff  et  ah  (2002)  only  selected  those  time 
delays  with  CC  values  larger  than  0.7  and  mean  coherences  above  0.7.  The  selection  of  an  "optimum"  threshold 
value  is  important  but  difficult.  If  it  is  set  too  high,  then  only  a  limited  number  of  very  accurate  differential  time  data 
are  available  to  constrain  the  relative  positions  of  earthquakes.  On  the  other  hand,  if  the  threshold  value  is  set  too 
low,  then  many  unreliable  differential  time  estimates  are  used  and  will  negatively  affect  the  final  relocation  results. 
Also  the  time  delay  estimate  calculated  with  the  CC  technique  may  not  be  reliable  when  the  noise  sources  in  the  two 
time  series  are  correlated. 

The  WCC  approach  and  its  variants  work  in  the  second-order  spectral  domain.  When  the  underlying  signal  inside 
two  similar  time  series  can  be  regarded  as  a  non-Gaussian  process  and  the  noise  sources  as  zero-mean  Gaussian 
processes,  the  similarities  between  the  two  time  series  can  be  better  compared  in  the  third-order  or  bispectrum 
domain  (Nikias  and  Raghuveer,  1987;  Nikias  and  Pan,  1988;  Yung  and  Ikelle,  1997).  This  bispectrum  method  takes 
advantage  of  the  fact  that  for  Gaussian  processes  only  all  spectra  of  order  higher  than  two  are  identically  zero.  Thus 
the  effect  of  correlated  Gaussian  noise  is  completely  suppressed  when  we  compare  the  two  time  series  in  the 
bispectrum  domain,  while  it  may  make  WCC  techniques  fail  to  work  well.  Noise  at  a  station  for  different  events  can 
be  expected  to  be  partially  correlated  due  to  a  combination  of  constant  predominant  noise  sources  with  time-varying 
amplitude  (microseisms,  wind  or  cultural  noise)  and  site  response  effects.  Cycle  skips  are  another  potential  WCC 
pitfall  that  bispectrum  analysis  may  help  detect. 

Figure  1  demonstrates  the  application  of  the  CC  and  bispectrum  methods  to  the  waveforms  of  two  closely  spaced 
M2.5  New  Zealand  earthquakes.  The  waveforms  for  these  events  recorded  at  nearby  stations  are  highly  similar,  but 
the  signal  to  noise  ratios  (SNRs)  for  the  records  at  station  MOW  (epicentral  distance  of  93  km)  are  noticeably  lower 
(Figs.  1 A  and  IB).  The  two  bispectrum-correlation  series  (Figs.  IE  and  IF)  peak  at  lag  6  and  7  respectively,  whereas 
the  two  CC  series  (Figs.  1C  and  ID)  reach  maxima  at  lags  very  different  from  each  other  and  also  from  the  two 
bispectrum-correlation  series.  Careful  examination  reveals  that  the  two  CC  series  reach  local  maxima  at  lag  7.  They 
fail  to  become  global  maxima  there,  however,  apparently  due  to  the  contamination  of  correlated  noise.  Figures  1G 
and  1H  plot  the  stacked  signals  for  both  raw  and  band-pass  filtered  waveforms.  In  both  cases  the  lower  trace  (the 
stacked  signal  obtained  by  shifting  the  waveform  of  the  first  event  with  bispectrum-determined  lags)  has  a  larger 
root-mean-square  (RMS)  amplitude  than  the  upper  trace,  which  is  obtained  after  shifting  the  waveform  of  the  first 
event  with  the  CC-determined  lags.  The  above  example  corroborates  the  findings  of  Nikias  and  Pan  (1988)  and 
Yung  and  Ikelle  (1997)  that  the  bispectrum  method  works  better  than  CC  techniques  in  getting  relative  time  delay 
between  two  similar  signals  contaminated  by  correlated  noise. 

In  our  approach,  we  compute  two  more  time  delay  estimates  between  the  waveform  pairs  with  the  bispectrum 
method  (the  raw  data  and  band-pass  filtered  data)  and  use  them  to  verify  the  WCC-determined  one  using  the  band¬ 
pass  filtered  data.  This  technique  can  provide  quality  control  over  the  selected  time  delay  estimates  and  potentially 
provide  more  differential  travel  time  data  for  close  events  pairs,  by  verifying  the  reliability  of  differential  times  that 
might  not  meet  the  threshold  criterion. 

The  process  of  time  delay  calculation  and  bispectrum  verification  can  be  briefly  described  as  follows.  For  the  k-th 
mutual  station  that  provides  waveforms  for  two  earthquakes  i  and  j,  we  rely  on  the  catalog  phase  (P  or  S)  picks  to 


form  the  data  windows  used  for  time  delay  calculations.  Then  using  the  band-pass  filtered  waveforms,  we  calculate 
both  the  CC  value  CQjk  and  the  relative  time  delay  in  lag  number  Ajjk(cc)  for  the  event  pair.  If  the  value  of  CQjk  is 
larger  than  a  threshold  CCsub,  we  also  perform  sub-sample  time  delay  calculation  to  get  Ayk(sub)  by  a  weighted  linear 
fitting  of  the  cross  spectrum  phase  (Poupinet  et  al.,  1984).  After  we  make  the  CC  calculations  for  all  the  mutual 
stations  between  events  i  and  j,  we  can  obtain  the  maximum  CC  value  CCymax  for  the  event  pair.  If  CCjjmax  is  larger 
than  a  second  threshold  CCbs,  we  will  perform  bispectrum  time  delay  estimation  for  each  mutual  station  of  the  event 
pair.  We  obtain  two  estimates  of  time  delay  in  lag  number,  i.e.  Aijk(bsl)  using  the  band-pass  filtered  waveforms  and 
Aijk(bs2)  with  the  raw  waveforms.  Thus  depending  on  the  extent  of  waveform  similarity  for  an  event  pair,  which  is 
measured  by  the  size  of  CQjmax,  four  possible  time  delay  estimations  are  carried  out  at  a  mutual  station  for  either  P 
or  S  waves. 

Besides  the  aforementioned  threshold  valud  that  researchers  often  use  to  select  time  delay  estimates,  which  we  term 
CChml,  three  other  parameters  CClim2,  CClim3  and  AUm  control  the  bispectrum  verification  process.  Generally  CChm2 
a  CChml  &  CClun3.  After  we  make  the  CC  calculations  across  all  the  mutual  stations  of  an  event  pair,  we  compare  the 
maximum  CC  value  CCymax  with  CChm2.  If  CQjmax  is  larger  than  CClim2,  we  will  verify  the  CC  time  delay  estimates 
for  those  stations  with  CCjjk  larger  than  CChm3.  A  CC  time  delay  estimate  Ayk(cc)  is  trusted  (or  passes  the  verification) 
if  its  differences  from  both  Ayk(bsl)  and  Ayk(bs2),  the  ones  determined  with  the  bispectrum  method,  do  not  exceed  the 
tolerance  limit  Ahm.  In  other  words,  for  an  event  pair  with  high  CQjmax  value  we  will  select  the  CC  time  delay 
estimates  as  long  as  they  pass  the  bispectrum  verification,  even  though  the  CC  values  associated  with  some  of  them 
are  smaller  than  CChml  and  would  not  be  used  under  the  threshold  criterion  often  adopted  by  other  researchers.  If 
CClunl  ^  CCijmax  <  CChm2,  we  will  make  the  verifications  only  for  those  stations  with  CCyk  larger  than  CChml.  If 
CCymax  falls  below  CChml,  the  CC  time  delays  for  the  event  pair  are  simply  discarded. 

As  an  example,  we  apply  this  technique  to  obtain  bispectrum-verified  WCC  differential  times  for  a  set  of  New 
Zealand  earthquakes,  and  relocate  521  of  the  events  using  the  DD  algorithm  hypoDD  (Waldhauser  and  Ellsworth, 
2000).  We  find  that  the  bispectrum-verified  time  delays  provide  more  clustered  earthquake  relocation  results  with 
lower  arrival  time  residuals  compared  to  the  threshold  criterion.  Figure  2  shows  the  relocation  results  for  a  subregion 
near  Lake  Wairarapa,  New  Zealand. 

Generally  for  an  earthquake,  fewer  S  arrivals  are  available  in  the  phase  catalog  than  P  arrivals  because  they  are  more 
difficult  to  select.  We  can  obtain  additional  bispectrum-verified  S  differential  times  for  those  waveforms  lacking 
catalog  S  picks,  using  the  predicted  S  arrivals  either  from  a  velocity  model  (Shearer,  1997)  or  from  the  S-P  time  of  a 
nearby  event.  Our  work  demonstrates  that  relocation  results  can  be  further  improved  with  these  additional  time 
constraints. 

Hierarchical  clustering  and  location 

We  have  implemented  and  tested  a  two-step  location  technique  that  offers  several  advantages  over  conventional 
techniques  due  to  its  use  of  WCC  and  DD  location  methods.  Step  1  event  location  is  based  solely  on  high-quality 
GT  event  waveforms  without  use  of  catalog  picks.  Step  2  location  estimates  offer  the  same  improvement  in  location 
accuracy  as  documented  by  the  use  of  high-resolution  DD  techniques  over  standard  techniques  (Waldhauser  and 
Ellsworth,  2000).  We  evaluated  the  technique  using  640  events  from  N.  California  (magnitude  range  4.3-S.3)  using  a 
6°x  6°  study  region  and  7  master  stations  from  the  UC  Berkeley  NCSN  array  (Fig.  3).  In  step  1,  event  locations  are 
estimated  via  CC  of  Hilbert-envelope  waveforms.  A  test  event  is  cross-correlated  against  a  suite  of  reference  events 
binned  by  latitude,  longitude,  and  depth.  The  step  1  location  derived  from  the  CC  global  maximum  is  then  used  as  a 
starting  location  in  the  step  2  refinement  using  the  DD  relocation  method  (hypoDD). 

All  events  were  initially  binned  using  the  catalog  locations  into  0.1°  latitude  and  longitude  bins  and  1 .5  km  depth 
bins.  We  selected  the  highest  magnitude  event  per  bin  as  a  reference  event.  We  assembled  a  suite  of  test  events  by 
selecting  the  second  highest  magnitude  event  from  each  multiple-event  bin.  After  data  quality  considerations,  this 
yielded  161  reference  events  and  19  test  events  for  use  in  the  step  1  waveform  location  technique  with  location 
accuracy  limited  to  the  binning  interval. 

Relocation  with  hypoDD  using  NCSN  data  from  25  Berkeley  network  stations  provided  GT  hypocenter  information 
for  the  test  events.  We  explored  4  different  parameter  settings  controlling  the  extent  of  clustering  and  residual  outlier 
rejection.  Average  standard  deviations  of  GT  locations  for  the  19  test  events  indicate  good  control  on  latitude  and 
longitude  (±0.015°  and  ±0.02°)  and  moderate  control  on  depth  (±1.23  km).  We  chose  GT  solutions  with  amean 
errors  (absolute  value)  of  0.03°,  0.04°,  and  1.26  km. 


In  the  step  1  procedure,  the  entire  seismogram  was  filtered  and  an  eigenvector  decomposition  of  the  covariance 
matrix  was  applied  to  the  3-component  data  (Rowe  et  al.,  2002)  prior  to  the  calculation  of  Hilbert-envelopes  for  the 
principal  direction.  Test  event  data  for  each  of  the  master  stations  were  cross-correlated  against  161  GT  events  with 
the  position  of  the  CC  global  maximum  determining  the  step  1  location  (Fig.  3).  The  Step  1  procedure  located  14  of 
19  test  events  to  within  20  km  (approximately  1  bin)  of  the  NCSN  catalog  location  with  location  error  for  the 
remaining  events  ranging  from  25-50  km.  Errors  in  test-event  depth  were  <  1.5  km  (1  depth  bin)  for  half  of  the 
events  and  on  the  order  of  2.5-10  km  for  the  remaining  step  1  relocated  events. 

The  step  2  procedure  attempts  to  improve  upon  the  step  1  epicentral  and  depth  estimates  based  on  the  use  of  catalog 
and  WCC  differential  times  with  hypoDD.  For  each  test  event,  subsets  of  the  catalog  and  WCC  data  are  prepared  for 
use  in  hypoDD  with  the  step  1  location  providing  the  initial  test-event  location. 

We  augmented  the  existing  P  and  S  picks  from  the  network  catalog  using  skew  and  kurtosis  analysis  followed  by  a 
manual  evaluation.  Many  catalog  P  picks  required  manual  adjustment.  We  then  prepared  WCC  data  based  on  3-  and 
5-second  P  windows  (Rowe  et  al.,  2002).  A  subset  of  non-test  event  data  (P  and  S  NCSN  phase  catalog  times  and 
WCC  P-pick  adjustments)  augmented  with  test  event  data  allows  clustering  of  the  test  event  with  adjacent  reference 
events.  Clustering  parameters  for  hypoDD  are  chosen  to  improve  vertical  control  by  strongly  weighting  pick  data 
from  event  pairs  with  larger  separation  distances  and  WCC  data  from  event  pairs  with  smaller  separation  distances 
(Waldhauser,  2001). 

We  evaluated  various  ways  of  constructing  pick  and  WCC  data  subsets.  Relocation  results  were  not  influenced  by 
the  choice  of  distance  cutoff  in  the  cluster  analysis.  We  used  a  cutoff  value  of  1°  (reference  event  to  test  event)  in 
our  construction  of  pick  and  WCC  data  subsets.  Construction  of  WCC  differential  time  data  subsets  also  depended 
on  quality  information,  absolute  arrival-time  difference  cutoff  (5-15  s)  and  the  minimum  CC  value  (0.7-0.9).  Both 
arrival-time  and  WCC  data  were  influenced  by  choice  of  clustering  parameters. 

We  also  experimented  with  use  of  bispectrum-correlation  (BC)  values  to  validate  WCC  information  with  low  CC 
values.  An  event  was  eligible  for  BC  analysis  if  a  station  CC  value  was  greater  than  a  threshold  value  cutoff.  Data 
from  a  subset  of  stations  per  event  were  selected  for  BC  calculation  based  on  a  minimum  CC  value  cutoff.  If  the  BC 
and  CC  lag  values  were  within  a  time-difference  tolerance  (0.25  s),  BC  data  was  allowed  into  the  DD  data  set.  The 
lowering  of  threshold  value  cutoffs  (0.5-0.7)  increased  inclusion  of  BC  data  from  10%  to  20%  and  led  to  a  slight 
improvement  in  step  2  depth  estimation  over  use  of  P  pick  WCC  data  with  standard  parameters  (0.7  CC  cutoff  and 
maximum  travel-time  difference  of  10  sec).  Although  we  have  not  at  this  point  adjusted  catalog  S  picks,  S-pick 
WCC  data  based  on  30-s  windowing  were  added  to  the  P-pick  WCC  data  as  a  final  effort  to  improve  depth 
estimation. 

To  test  the  stability  of  the  step  2  procedure,  locations  were  calculated  using  a  range  of  starting  depths  (1.0,  2.5,  5.0, 
10.0,  15.0,  and  20.0  km)  for  each  test  event.  Location  mean  and  standard  deviation  values  (Fig.  4)  show  that  depths 
are  less  constrained  for  events  outside  of  the  ring  of  master  stations,  in  particular  the  test  events  positioned  near 
34.3°.  This  is  indicated  by  the  tendency  of  the  mean  locations  to  shift  towards  the  center  of  the  test  depth  range  and 
by  larger  depth  standard  deviations.  Average  standard  deviation  values  of  latitude,  longitude,  and  depth  for  test 
events  near  34.3°  are  larger  (0.027°,  0.019°,  4.6  km)  than  within  station  coverage  (0.004°,  0.013°,  and  2.35  km) 
suggesting  tradeoff  in  latitude  and  depth  for  the  outside  events.  We  attempted  to  improve  the  stability  of  the  outlier 
events  by  adjusting  WCC  P-window  length,  including  BC  data,  and  adjusting  parameters  that  control  the  use  of 
WCC  and  BC  data.  A  5-second  P  window  with  a  0.9  cross-correlation  cutoff  was  used  to  produce  the  results  shown 
in  Figure  4.  The  use  of  the  BC  data  or  other  WCC  parameter  settings  yields  similar  results. 

We  included  S-wave  WCC  data  to  calculate  step  2  relocations  displayed  in  Figure  5,  as  it  slightly  lowered  relocation 
errors  for  some  of  the  events  outside  station  coverage  (34.3°  latitude).  The  hypoDD  relocation  improves  both 
epicenter  and  depth  estimates  (Fig.  5)  both  in  resolution  (GT  error  <  bin  intervals)  and  in  accuracy.  Seventeen  of  19 
test  events  have  hypocenter  location  errors  <  15  km,  with  1 1  events  relocated  to  within  5  km  of  GT.  Depth  errors  for 
16  events  are  within  2.5  km  of  GT.  Mean  latitude,  longitude  and  depth  errors  (absolute  value)  for  test  events  outside 
station  coverage  are  larger  (0.059°,  0.051°,  and  2.5  km  depth)  by  roughly  a  factor  of  2  than  mean  errors  within 
station  coverage  (0.028°,  0.031°,  and  1.02  km  depth).  Step  2  location  mean  errors  for  events  within  station  coverage 
are  within  the  mean  location  error  of  the  GT  solution.  Improvement  of  the  S  arrival-time  picks  over  catalog  times 
should  lead  to  improved  CC  data  and  depth  estimation.  The  accuracy  of  the  relocation  estimates  even  with  coarse 
master  station  coverage  suggests  this  approach  could  have  broader  application  outside  of  this  particular  GT  dataset 
study  region. 


Double-difference  tomography 


We  have  developed  a  new  DD  tomography  method  (Zhang  and  Thurber,  2003),  which  takes  a  step  further  of  the  DD 
location  method  (Waldhauser  and  Ellsworth,  2000)  to  solve  for  3D  velocity  structure  simultaneously  with 
hypocenter  locations  using  both  catalog  picks  and  differential  arrival  times  (WCC  and  catalog).  At  the  local  scale 
(10s  to  100s  of  kilometers),  the  earth  can  be  represented  with  a  flat  model.  We  use  the  pseudo-bending  ray-tracing 
algorithm  (Um  and  Thurber,  1987)  to  find  the  rays  and  calculate  the  travel  times  between  events  and  stations.  The 
model  is  represented  as  a  regular  set  of  3D  nodes  and  the  velocity  values  are  interpolated  by  using  the  tri-linear 
interpolation  method.  The  hypocentrer  partial  derivatives  are  calculated  from  the  direction  of  the  ray  and  the  local 
velocity  at  the  source.  The  ray  path  is  divided  into  a  set  of  segments  and  the  model  partial  derivatives  (calculated  in 
terms  of  fractional  slowness  perturbation,  so  that  the  derivatives  are  related  to  path  length)  are  evaluated  by 
apportioning  the  derivative  to  its  8  surrounding  nodes  according  to  their  interpolation  weights  on  the  segment 
midpoint  (Thurber,  1983). 

We  have  tested  the  local-scale  version  of  DD  tomography  (tomoDD)  with  a  synthetic  dataset  that  was  constructed 
on  the  basis  of  an  idealized  model  of  the  velocity  structure  of  the  San  Andreas  Fault  in  central  California  (Kissling  et 
al.,  1994)  shown  in  Figure  6.  The  events  and  stations  used  to  construct  the  synthetic  dataset  are  from  the  actual 
seismicity  and  USGS  stations  in  the  Loma  Prieta  region  (Fig.  7).  We  added  Gaussian  random  noise  with  zero  mean 
and  a  standard  deviation  of  0.04  s  to  all  the  true  arrival  times.  In  addition,  we  also  added  a  constant  noise  term  to  the 
arrivals  at  each  station  from  a  uniform  distribution  between  -0.3  s  and  0.3  s.  This  simulates  the  case  that  the 
systematic  errors  (model  errors  and  pick  bias)  associated  with  the  arrival  times  are  larger  than  the  random  ones.  We 
construct  the  pseudo  ’'cross-correlation"  data  directly  from  the  absolute  arrival  times  by  differencing  the  synthetic 
arrival  times  at  common  stations  for  events  pairs  within  20  km.  As  a  result,  the  "cross-correlation"  data  are  more 
accurate  than  the  absolute  data. 

We  first  use  the  DD  location  algorithm  hypoDD  to  relocate  the  events  (Waldhauser,  2001).  The  ID  velocity  model 
used  is  the  Dietz  and  Ellsworth  (1990)  model,  based  on  seismic  refraction  and  earthquake  data.  We  use  both  the 
residual  weighting  and  distance  weighting  schemes  in  the  inversion  to  control  the  large  residuals,  as  did  Waldhauser 
and  Ellsworth  (2000).  The  absolute  differences  between  this  set  of  event  relocations  and  the  true  locations  and  their 
standard  deviations  are  shown  in  Table  1.  The  results  indicate  that  the  event  relocations  from  the  DD  location 
algorithm  based  on  a  ID  velocity  model  have  a  substantial  bias  (>1  km  in  each  coordinate  direction)  from  the  true 
locations.  This  bias  is  caused  by  a  difference  between  the  velocity  model  (horizontal  layers)  used  to  calculate  the 
DD  locations  and  the  true  velocity  model  (vertical  "sandwich")  used  to  generate  the  data.  The  heterogeneity  of  the 
true  velocity  model  makes  path  anomalies  different  for  different  events.  However,  the  DD  location  method  assumes 
that  the  path  anomalies  are  location-independent,  and  this  assumption  introduces  bias  into  event  locations  (Wolfe, 
2002). 

Next,  we  show  the  results  of  applying  the  standard  tomography  method  to  the  relatively  noisy  absolute  arrival  times. 
The  inversion  starts  from  the  same  ID  velocity  model  as  the  DD  location  method.  The  X-Y  nodes  used  to  represent 
the  velocity  model  are  shown  in  Figure  7;  in  depth,  nodes  are  positioned  at  0,  3,  7,  1 1,  and  16  km.  The 
computational  algorithm  is  identical  to  that  for  DD  tomography,  but  the  differential  times  are  excluded.  Both 
damping  and  smoothing  (with  weight  5.0  in  both  the  horizontal  and  vertical  directions)  are  applied  to  the  inversion 
to  make  the  solution  more  stable.  Figure  8a  shows  horizontal  slices  through  the  velocity  structure  obtained  from  the 
standard  tomography.  The  absolute  difference  between  this  velocity  model  and  the  true  velocity  model  has  a  median 
value  of  0.164  km/s,  a  mean  value  of  0.245  km/s  and  a  standard  deviation  of  0.249  km/s.  The  main  features  of  the 
true  velocity  model  are  evident  in  the  depth  slices.  The  event  locations  have  smaller  errors  than  those  from  the  DD 
location  method  (Table  1). 

DD  tomography  uses  both  the  noisy  absolute  arrival  times  and  more  accurate  differential  times.  We  also  apply 
damping  and  the  smoothing  constraints  to  stabilize  the  solution,  with  the  residual  weighting  used  to  penalize  the 
large  residuals  during  the  inversion.  DD  tomography  starts  with  the  same  ID  velocity  model  and  uses  the  same 
inversion  grid  as  the  standard  tomography.  Figure  8b  shows  the  horizontal  slices  of  the  velocity  structure  from  DD 
tomography.  The  absolute  difference  between  this  velocity  structure  and  the  true  velocity  structure  has  a  median 
value  of  0.136  km/s,  a  mean  value  of  0.178  km/s,  and  a  standard  deviation  of  0.164  km/s.  DD  tomography 
characterizes  well  the  low  velocity  zone  of  the  true  velocity  model,  especially  at  the  depths  of  3  km  and  7  km. 
Subtracting  the  DD  tomography  solution  and  the  standard  tomography  solution  from  the  true  model  (Fig.  8c  and  d), 
we  find  that  the  velocity  model  from  DD  tomography  has  a  more  correct  value  in  the  low-velocity  trough  except  at 
the  depth  of  0  km,  where  the  ray  paths  from  event  pairs  almost  completely  overlap  and  the  accuracy  of  velocity 


structure  is  mainly  controlled  by  absolute  catalog  data.  This  indicates  that  DD  tomography  recovers  better  the  low- 
velocity  zone,  and  overall  the  velocity  model  is  more  accurate  than  the  standard  tomography.  Compared  with  the 
standard  tomography  method,  DD  tomography  also  produces  more  accurate  event  locations  (Table  1). 

At  the  regional  scale  (100’s  to  1000's  of  km),  however,  the  earth  should  not  be  treated  as  flat,  but  as  a  sphere.  Some 
velocity  discontinuities  such  as  Conrad,  Moho,  and  subducting  slab  boundary  should  also  be  taken  into  account. 
Finite-difference  ray  tracing  algorithms  developed  by  both  Podvin  and  Lecomte  (1991)  and  Hole  and  Zelt  (1995)  are 
able  to  accurately  calculate  first-arrival  times  in  the  presence  of  extremely  severe,  arbitrarily  shaped  velocity 
discontinuities.  These  algorithms  solve  a  finite  difference  approximation  to  the  Eikonal  equation  in  the  regularly 
gridded  velocity  structure  through  a  systematic  application  of  Huygens'  principle.  This  procedure  can  explicitly  take 
into  account  different  propagation  modes  including  transmitted  and  diffracted  body  waves,  and  head  waves  in 
addition  to  direct  waves.  The  principle  of  reciprocity  is  utilized  by  placing  the  source  point  at  the  seismic  station 
location  and  interpolating  the  travel  times  computed  at  the  grid  nodes  to  match  the  exact  earthquake  location 
(Flanagan  et  al.,  1999). 

To  adapt  the  finite-difference  algorithm  to  the  regional  scale,  the  curvature  of  the  Earth  must  be  taken  into  account. 
Earth  flattening  transformation  is  only  valid  for  ID  velocity  models  and  not  for  3D  velocity  models.  Following 
Flanagan  et  al.  (2000),  we  solved  this  problem  by  parameterizing  a  spherical  surface  inside  the  Cartesian  volume  of 
grid  nodes.  The  coordinate  center  was  put  at  the  surface  of  the  Earth,  positive  X  and  Y-axes  point  to  the  direction  of 
North  and  West,  and  positive  Z-axis  points  downward,  respectively.  The  grid  nodes  above  the  Earth  surface  (air 
nodes)  are  attributed  to  the  true  velocity  values  traveling  in  the  air.  As  a  result,  all  the  rays  travel  inside  the  Earth. 
The  regular  uniform  velocity  (slowness)  grid  nodes  used  to  calculate  the  travel  time  field  using  finite-difference 
algorithms  are  interpolated  from  the  non-uniform  inversion  grid  nodes  through  tri-linear  interpolation.  If  an 
inversion  grid  node  is  2  km  above  the  Earth's  surface,  then  it  is  treated  as  an  air  node  and  its  value  is  fixed  during  the 
inversion.  We  treat  each  station  as  a  source  and  calculate  travel  times  to  all  velocity  nodes  in  the  volume.  The  travel 
time  from  a  station  to  each  earthquake  is  interpolated  from  its  8  neighboring  nodes  through  tri-linear  interpolation. 
The  ray  path  from  the  earthquake  to  the  station  is  found  iteratively  with  increments  opposite  to  the  travel  time 
gradient.  After  these  setups,  DD  tomography  is  adapted  to  the  regional  scale  with  ray  paths  and  related  partial 
derivatives  calculated  by  the  finite-difference  method.  We  have  tested  the  regional-scale  version  of  DD  tomography 
on  data  from  the  subduction  zone  beneath  northern  Honshu,  Japan,  and  obtained  the  velocity  structure  with 
unprecedented  resolution. 

CONCLUSIONS  AND  RECOMMENDATIONS 

(1)  WCC  is  vulnerable  to  problems  associated  with  correlated  noise  that  can  be  overcome  using  bispectrum  analysis. 
It  is  our  recommendation  that  a  verification  approach  comparable  to  that  described  here  be  adopted  for  earthquake 
location  applications  for  whiph  very  high  accuracy  is  essential. 

(2)  A  hierarchical  location  approach,  for  "new"  events  relative  to  a  set  of  reference  events,  that  starts  with  CC  of 
long  records  at  low  frequencies  to  identify  the  approximate  source  region  based  solely  on  waveform  similarity  and 
then  switches  to  DD  location  using  CC  of  short  records  at  high  frequencies  works  effectively  for  a  regional  dataset 
from  California.  This  approach  should  next  be  tested  in  a  real-time  monitoring  environment. 

(3)  Algorithms  for  DD  tomography  have  been  successfully  developed  for  both  local  (10's  to  a  few  100  km)  and 
regional  (100’s  to  1000’s  of  km)  scale.  Synthetic  tests  show  the  method  yields  locations  with  greater  accuracy  than 
either  DD  location  or  standard  tomography,  and  a  velocity  model  that  is  closer  to  the  true  model  than  standard 
tomography.  The  DD  tomography  approach  works  with  either  waveform-based  or  catalog-based  differential  time 
data,  and  thus  is  applicable  to  any  dataset  to  which  standard  tomography  can  be  applied. 
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Figure  1.  (A)  Windowed  waveforms  (1  second  before  and  2.82  seconds  after  the  catalog  P  pick)  for  two  New 
Zealand  earthquakes  recorded  at  station  MOW.  The  station  sampling  rate  is  50  sample/sec  and  the  signals  are 
aligned  by  the  P  picks.  (B)  Windowed  waveforms  band-pass  filtered  between  1  and  6  Hz.  (C,  D)  CC  series  for  the 
(C)  unfiltered  (raw)  and  (D)  filtered  waveforms.  (E,F)  Bispectrum-correlation  series  for  the(E)  raw  and  (F)  filtered 
waveforms.  (G)  Stacked  raw  waveforms.  The  upper  trace  is  obtained  by  shifting  the  waveform  of  the  first  event  with 
the  CC  determined  lag  relative  to  the  second  one,  while  the  lower  trace  is  computed  after  shifting  the  waveform  of 
the  first  event  with  the  lag  calculated  with  the  bispectrum  method.  The  root-mean-square  amplitudes  for  the  two 
traces  are  0.25  and  0.27  respectively.  (H)  Same  as  G  for  the  stacked  band-pass  filtered  waveforms.  The  root-mean- 
square  amplitudes  for  the  two  traces  are  0.32  and  0.34  respectively. 


Figure  2.  Locations  of  53 
earthquakes  near  lake 
Wairarapa,  NZ.  (A)  Before 
relocation;  (B)  relocated  with 
catalog  differential  travel 
times;  (C)  relocated  with  CC 
differential  travel  times 
chosen  with  the  threshold 
criterion;  (D)  relocated  with 
CC  differential  travel  times 
verified  with  the  bispectrum 
method;  (E)  relocated  with 
both  catalog  and  threshold- 
chosen  CC  differential  travel 
times;  (F)  relocated  with  both 
catalog  and  bispectrum- 
verified  CC  differential  travel 
times. 


Figure  3.  A)  Map  of  study  region  showing  the  7  master  stations  used  in  the  step  1  procedure,  161  reference  events, 
and  19  test  events.  Stations  are  indicated  by  triangles  and  labeled  by  their  BK  network  station  abbreviation.  Test  and 
reference  events  are  indicated  by  **’  and  *+’  symbols  respectively.  Usage  of  hypoDD  to  construct  GT  for  test  events 
and  in  the  step  2  relocation  includes  catalog  data  from  18  Berkeley  network  stations  in  addition  to  the  7  master 
stations  for  a  total  station  range  of  .35.9°  to  41.9°  latitude  and-1 19.3°  to  -124.07°  longitude.  B)  Step  1  derived 
mislocation  of  hypocenter  indicated  by  lines  connecting  the  step  1  location  with  the  hypoDD  based  GT  location  ('o' 
symbol).  Events  positioned  outside  master  station  coverage  (near  34°  latitude)  have  larger  mislocation  errors  relative 
to  GT  than  test  events  within  the  circumference  of  master  station  coverage. 


Figure  4.  Stability  test  showing 
poorer  depth  constraint  for  events 
outside  the  set  of  master  stations. 
The  mislocation  is  indicated  by 
solid  line  connecting  the  mean  step 
2  location  with  the  hypoDD  GT 
location  (’o’  symbol).  The  vertical 
and  horizontal  lines  at  the  mean 
step  2  location  indicate  standard 
deviations  of  the  step  2  estimates 
of  depth  and  latitude. 


Figure  5.  Step  2  mislocation  of 
hypocenters  indicated  by  lines 
connecting  the  step  2  location 
with  the  hypoDD  based  GT 
location  ('o'  symbol). 
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Figure  6.  A  horizontal  slice  through  the 
true  synthetic  velocity  model.  The  true 
velocity  model  in  3D  is  similar  to  a 
"vertical  sandwich,"  with  velocity 
constant  with  depth. 


Figure  7.  Event  locations  (filled  circles)  and  stations  (triangles) 
used  for  the  synthetic  data  set.  The  inversion  grid  used  in  the 
standard  and  DD  tomography  solutions  is  shown  as  crosses.  The 
inversion  grid  points  are  at  X=-35,  -15,  0,  2, 4,  6,  20,  35  km,  at 
Y=-60,  -40  -20,  0,  20,  40  km,  and  at  Z=0,  3,  7,  11,  16  km. 
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Figure  8.  Horizontal  slices  through  the  velocity  models  from  (a)  standard  tomography  and  (b)  DD  tomography,  and 
the  difference  between  (c)  the  DD  tomography  solution  and  the  true  model  and  (d)  the  standard  tomography  solution 
and  the  true  model.  Black  dots  indicate  the  earthquake  hypocenters  within  half  die  grid-size  of  the  slice. 

Table  1.  The  absolute  differences  between  the  true  locations  and  those  from  the  DD  location  method  based  on  ID 
velocity  model,  standard  tomography,  and  DD  tomography. 


Median  value  (km) 


Standard  deviation  (km) 


